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Abstract. The interactions between the components of complex networks 
are often directed. Proper modeling of such systems frequently requires the 
construction of ensembles of digraphs with a given sequence of in- and out-degrees. 
As the number of simple labeled graphs with a given degree sequence is typically 
very large even for short sequences, sampling methods are needed for statistical 
studies. Currently, there are two main classes of methods that generate samples. 
One of the existing methods first generates a restricted class of graphs, then uses 
a Markov Chain Monte-Carlo algorithm based on edge swaps to generate other 
realizations. As the mixing time of this process is still unknown, the independence 
of the samples is not well controlled. The other class of methods is based on the 
Configuration Model that may lead to unacceptably many sample rejections due 
to self-loops and multiple edges. Here we present an algorithm that can directly 
construct all possible realizations of a given bi-degree sequence by simple digraphs. 
Our method is rejection free, guarantees the independence of the constructed 
samples, and provides their weight. The weights can then be used to compute 
statistical averages of network observables as if they were obtained from uniformly 
distributed sampling, or from any other chosen distribution. 



PACS numbers: 02.10.Ox, 02.50.Ey, 89.75.Hc, 07.05.Tp 
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1. Introduction and definitions 

In network modeling problems [XJ [2l El SJ [5j EJ [T] , one often needs to generate ensembles 
of graphs obeying a given constraint. A typical constraint is the case when the only 
information available is the degrees of the nodes, and not the actual connectivity 
matrix. Note that the node degrees by themselves, that is the degree sequence in 
general does not determine a graph uniquely: there can be a very large number 
of graphs having the same degree sequence [8]. Full graph connectivity is uniquely 
determined by the degree sequence only for a special class of sequences (see Ref. [S] 
for the case of undirected graphs). 

Often, the interest lies in the study of network observables, as determined by 
the given sequence of degrees, and unbiased by anything else. These can be graph 
theoretical measures, or properties of processes happening on the network (e.g., 
spreading processes, such as of opinion or disease). The problem of creating and 
sampling graphs with a given degree sequence, i.e., degree-based graph construction |101 
1 lj , is a well-known and challenging problem that has attracted considerable interest 
amongst researchers [U EQ] Ell El El El 113 El ED EEE1 OS HU HH1 Ull 
28, 29 , 30J. There are two main classes of algorithms that are used today to achieve the 
construction of graphs with given degree sequences. One of them is typically referred 
to as "switching" or edge-swap based (13] EH ESI 12Z| , while the other one is usually 
called "matching" or stub-matching based [1E1E1E3EI1EI1EI1I2E1I3I]- Switching 
methods repeatedly swap the ends of two randomly chosen edges within a Markov 
Chain Monte-Carlo (MCMC) scheme until a new, quasi-independent, sample is 
produced. Unfortunately, the mixing time of MCMC schemes for arbitrary sequences is 
not known in the general case. The other class consists of direct construction methods, 
which perform pairwise matchings of the half-edges emanating from randomly chosen 
nodes until all edges are realized. Unfortunately, this method can easily generate 
multiple edges and self-loops, i.e., edges starting and ending on the same node, after 
which the sample must be rejected in order to avoid biases [35] . For a comparison of 
the two classes of methods see Ref. |23| . 

Recently, a novel degree-based construction [TO] and sampling method [TT] was 
introduced for undirected graphs, which has a worst-case scaling of O(NM), where 
M is the number of edges (2M is the sum of the degrees, which are given). A similar 
method was obtained independently in Ref. [34] , but that method is less efficient, 
with a worst-case scaling of 0(N 2 M). Although the algorithm in Ref. is a 
direct construction method using stub-matchings, it is rejection free, the samples are 
statistically independent and the algorithm also provides a weight for every realization. 

In many systems the interaction between two entities is not mutual but has 
a direction from one to the other, such as in the cases of human relationships in 
social networks [36], gene interactions in regulatory networks, trophic interactions in 
food webs |20l 121] , etc. Such systems require a representation by directed graphs 
(digraphs). In fact, undirected graphs can be interpreted as digraphs in which there 
are two, oppositely directed edges for each connected pair of nodes. Here we present 
a generalization of the degree-based graph construction problem to directed graphs. 
Some of the necessary mathematical foundations, laid down in Ref. |30_, are here used 
and expanded to introduce a digraph construction and sampling algorithm. Although 
the approach follows closely the one introduced by us for the undirected case [TT], the 
generalization is not at all straightforward, and there are significant differences that 
the directed nature of the links induces. 
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Figure 1. Examples of realizations of graphical bi-degree sequences. Panels a) 
and b) show two non-isomorphic realizations of the same bds. Panel c) shows 
a digraph that cannot be obtained via the Havel-Hakimi algorithm for digraphs. 
Panel d) shows a realization of a different bds. Panel e) illustrates that not all 
possible connections lead to a simple digraph even if a bds is graphical: in fact, 
the connections in the figure break the graphical character. 



Before we present our algorithm, we introduce some notations, based on Ref [30 . 
Let us denote by and the in- and out-degrees of a node i. Given the sequence 
D = | (d^, d^^j , (d^, d^^j , . . . , (d$ , d^^j j of non-negative integer pairs, we want 

to construct a simple directed graph G(V, E) such that node k £ V has (d^ \d^) 
for its in- and out-degrees, respectively, for all k = 1, 2, . . . , N. A simple directed 
graph is a graph that has no self-loops, nor multiple directed edges in the same 
direction between two nodes. There can be at most two edges between a pair of 
nodes, oppositely directed. We call the sequence D a bi-degree sequence (bds for 
short). When there is a simple digraph with a given bds D for its degrees, we say that 
the bds is graphical and that the digraph realizes D. Equivalcntly, we will also talk 
about "graphically" as a property. We distinguish realizations as labeled digraphs, 
and do not deal here with isomorphism questions. That is, if two realizations are 
identical up to a permutation of their indices, i.e., they are isomorphic, we will still 
consider them distinctly. In order to avoid isolated nodes, in the following we will 
assume that + d^ > 0, for all j = 1, . . . , N. As examples, Figs. Fk) and IJd) show 
two realizations of the bds Di = {(1,0), (1,2), (2,2), (2, 1), (0,1)}, and Fig. [if) shows 
a realization of D 2 = {(3, 0), (3, 0), (1, 2), (1, 2), (1, 2), (1, 2), (1, 2), (1, 2)}. Examples 
of non-graphical bds are the sequences D 3 = {(2, 2), (2, 1), (1, 3), (1, 1)} and D 4 = 
{(5, 6), (5, 6), (5, 6), (4, 3), (3, 3), (2, 1), (2, 1), (1, 1)}. 

Notice that even if a bds is graphical, not all connection sequences are guaranteed 
to end up with a simple digraph. For example, Fig. [TJi) shows a simple digraph 
realization of D5 = {(0, 1), (2, 0), (1, 2), (2, 2)}. However, if we were to place the first 
four edges as in Fig. [l£), we would break graphicality: from there on, we would not 
be able to complete the realization of the bds without creating either self-loops or 
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Figure 2. The construction of a digraph realizing a given bds proceeds by 
connecting the out-stubs of the nodes to the in-stubs of other nodes. In this 
"bipartite" representation the vertical grey bars represents single nodes. 



multiple edges. Hence, it is important to find an algorithm that builds digraphs with 
a given bds. As we will see, this is a challenging problem in itself. 

An algorithm that builds a digraph from a given bds sequentially connects the 
out-links of a node to the in-links of others. We can think of these out- and in-links 
as "out-stubs" and "in-stubs" emanating from a node, that are paired up with the 
corresponding stubs of other nodes. An intuitive representation of this is shown in 
Fig. [2] As the graph construction algorithm proceeds, the number of stubs of the nodes 
decreases. At any time during this process we will call the number of remaining in- 
stubs and out-stubs of a node its residual in- and out-degrees, and the corresponding 

bi-degree sequence D = | ^<ij l \ d'f^ , ( d^\ ■ ■ • , ( d$ , d$^j j the residual bds. 

Finally, another concept we will need to use in what follows is the notion of normal 
order [30 , which is essentially the lexicographic order on the bds. That is, we say that 

(i) (i) (i) 

a bds is in normal order, if for all 1 < j < N — 1, we have either dj > or, if dj — 

df +1 , then df > 4+i- Thus > the bds °6 = {(5,2), (4,4), (4,3), (2,5), (2,4), (2,1)} 
shown in Fig. [2j is arranged in normal order. Once a bds is in normal order, we will 
use the words 'left' or 'right' to describe the directions towards lower or higher index 
values in the sequence. 

The remainder of this paper is organized as follows: Section [2] introduces the 
fundamental mathematical notions and algorithmic considerations that are at the 
basis of our digraph construction algorithm. Section [3] presents the algorithm and 
its derivation details. Readers interested only in the algorithm itself may skip 
Subsection |3.1| and proceed to the summary described in the beginning of Section [3] 
and in Subsection 3.2 Section ^ deals in detail with the digraph sampling problem, 
provides the derivation of the sample weights and presents a simple example. Section [5] 
is dedicated to the complexity of the algorithm, and Section [6] concludes the paper. 



2. Mathematical foundations 



As seen from the examples above, not all sequences of non-negative integer pairs 
can be realized by simple digraphs. The sufficient and necessary conditions for the 
realizability of a bds are given by the "FR" theorem [371 1551 155] : 

Theorem 1 (Fulkerson-Ryser) A sequence of non-negative integer pairs D = 
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| (d[ i} , 4°') , . . . , f d$ , d ( °^ | with d[ i} >df>...> d$ is graphical iff 

df <N-1, d { ° ] < N - 1, 1 < i < N (1) 

N N 

£4° = £<< 0) . ( 2 ) 

i=l i=l 

and /or a// 1 < fc < N - 1 : 



k k N 



< E min { fc - Mi o) }+E min {M 0) } • (3) 

z— 1 i—1 i—k+1 

Given a bds, we can easily test if it is graphical using this theorem, and thus we will 
also refer to it as the "FR test". Condition Q states that both the number of in- 
and out-degrees for all nodes must be no larger than the number of other nodes it 
could connect to, or receive connections from. Condition Q is a consequence of the 
requirement that every out-stub must join an in-stub somewhere else; the sequence 
D3 given in one of the above examples is not graphical because it fails this condition. 
Condition (|3| is less intuitive. Its left hand side is the total number of in-stubs that 
the group of k highest in-degree nodes can receive. Within this group, a node's out- 
stubs can absorb no more of those in-stubs from the same group than its out-degree or 
k — 1 (it cannot absorb from itself), whichever is smaller (giving the first sum on the 
rhs of ([3j| ) . Outside of this group, a node cannot absorb more of those in-stubs than 
its out-degree or k, whichever is smaller (the second sum on the rhs of Q). Hence, 
the necessity of ^ . For the complete proof see Refs. [381 I3U] . Note that the example 
sequence D4 above fails condition [3] for k = 3. The FR test is the directed version of 
the Erdos-Gallai (EG) theorem (test) for undirected graphs. 

An important note is that bi-degree sequences are less constraining than 
undirected ones. The out-stub of a node is always connected to an in-stub of another, 
not affecting that node's out-stubs, whereas such distinction does not exist for the 
undirected case. Alternatively, if we disregard for a moment the directionality of the 
links and consider the degree of the node to be the sum of its in- and out-degrees, 
then the corresponding graph realizing the bds can have two edges running between 
the same pair of nodes, whereas this is not allowed in the undirected case. 



2.1. Algorithmic considerations 

The FR theorem only tests for graphicality, but it does not provide an algorithm for 
constructing the digraph(s) realizing the given bds. At first sight this might not seem 
an issue. However, the sequence D5 in Figs. [TJl and [IJ) reminds us that graphicality 
can easily be broken by a careless connection of stubs. Clearly, for the purposes of 
digraph construction, it should not matter which edges we create first, as long as we 
make sure that every connection made does not break graphicality. In other words, the 
possibility to create the rest of the edges, so that a simple digraph results in the end, 
must always be preserved. Thus, the key for the creation of an algorithm that builds 
simple digraphs realizing a given bds without rejections is in a theorem that allows us 
to check if we would break graphicality by placing a specific connection. Indeed, such 
theorems exist, and they will be discussed below. However, interestingly, they require 
that connections be made from the same node, until all its stubs are used away into 
edges. That is, assuming that we already made some connections from a given node 
i, preserving graphicality, these theorems give necessary and sufficient conditions for 
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keeping graphicality by the next connection still involving node i. Simply put, they 
won't work in general, if we would attempt a new connection from j to k, where 
j, k 7^ i, while node i still has dangling stubs. 

The connections already made from i to some set of nodes X^ represent a constraint 
for the new connections from i, as these novel connections must avoid the set X\. We 
call such a constraint associated to a node a star constraint on that node. Once all 
the stubs of node i are connected into edges while preserving graphicality, we obtain a 
graphical residual sequence D' on at most N — 1 nodes. Clearly, the new connections 
we make from this point on will not be constrained in any way by the connections we 
made from node i. For the purposes of realizing the sequence D' we can just simply 
remove node i with its fully completed connections, create a realization by a simple 
graph of D', then, in the end, add back node i with its connections to this graph in 
order to obtain a realization of D. The comments above hold both for the undirected 
and directed cases. 

One might think of using the EG test for the undirected case and the FR test for 
the directed case on a residual degree sequence to decide if graphicality was broken 
after attempting a new connection from the same node. For the undirected case, we 
have shown in Ref. [TT] that the passing of the EG test by the residual sequence is only 
a necessary condition, if there is already a star constraint on a node. For example, 
consider the graphical degree sequence d = {6, 5, 5, 3, 3, 2, 1,1}, and assume that we 
made connections from node i = 3 to nodes A3 = {1,6,7}. The residual sequence after 
these connections is d' = {5, 5, 2, 3, 3, 1, 0, 1}. It is easy to check that it passes the EG 
test. However, we will break graphicality with every realization of d', because it will 
form a double edge with one of the existing connections from node i = 3 to X s . Thus, 
additional considerations have to be made to ensure the graphicality of the residual 
sequence for the undirected described in [TT]. For the directed case here we 

use the sufficient and necessary conditions for graphicality under star constraints as 
provided by Theorem [2] below, proven in Ref. (30] . 

From now on, we will always talk about algorithms that first finish all the out- 
stubs of a node before moving onto another node with non-zero out-degree. In the 
case of a graphical bds, once all the out-degrees of all the nodes have been connected 
into directed edges, we are guaranteed to have completed a digraph, because the total 
number of in-stubs equals the total number of out-stubs, according to property Q . 

2.2. Theorems on which the algorithm is based 

An algorithm that builds graphical realizations of degree sequences of simple undirected 
graphs is the Havel-Hakimi (HH) algorithm [TUl ITT]: we choose any node with non- 
zero residual degree, then we connect all its stubs to nodes with the largest residual 
degrees avoiding self and multiple connections. This process is repeated with other 
nodes until all stubs of all nodes are used. There is a corresponding version of 
the HH algorithm for bi-degree sequences as well, introduced first in Ref. [T2], then 
rediscovered independently in Ref. [3U] , the latter providing an alternative proof. The 
HH algorithm for bds proceeds as follows: given a normal-ordered bds, choose any 
node with non-zero residual out-degree, then connect all its out-stubs to nodes with 
the largest residual in-degrees, without creating multiple edges running in the same 
direction, nor self-loops. Reorder in normal order the residual sequence and repeat 
this process until all stubs of all nodes are used. While for any given bds, the HH 
algorithm will construct a set of digraphs, it cannot construct all possible digraphs 
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realizing the same sequence, as shown in Ref. [30 . For example, the HH algorithm can 
never result in the digraph shown in Fig.JTj;) realizing the example sequence D 2 above. 
It is easy to see why: there are two kinds of nodes in this example, with bi-degrees 
(3, 0) and (1, 2). The only nodes with non-zero out-degrees are the (1, 2) types. Using 
the HH algorithm, we would have to connect both out-stubs of such a node to the 
nodes with the largest in-degrees, that is to the two (3, 0) types. However, the digraph 
in Fig. [TJi) does not have a (1,2) node being connected to both (3,0) nodes, yet it 
realizes the sequence. The limitation of the HH algorithm comes from the fact that 
it prescribes to connect the out-stub of a node i to an in-stub of the node with the 
largest residual in-degree that does not yet receive a connection from node i. However, 
there can be other nodes whose in-stubs can form a connection with an out-stub of i 
without breaking graphically. This shows the importance of finding a method able 
to build not just a realization of a bds, but all the possible realizations of any given 
bds. _ _ 

In the remainder, given a residual bds D, we denote by Ai (D) the allowed set of 
i, i.e., the set of all nodes to which an out-stub of i can be connected without breaking 
graphicality. Also, let us denote by Xi (D) the set of nodes to which connections were 
already made from i, thus representing the star constraint at that stage. 

The graphicality test under a star constraint on node i is provided as Theorem [2] 
below. In order to announce it, however, we need to introduce one more definition. 
Consider a bds D and a given node i with out-degree > from this bds. Let 

us also consider a subset of nodes S C V such that l^l < d!f\ where 15*1 denotes the 

(i) 

number of nodes in S, i.e., its size, and for every node j £ S, dj > 0. Next, we take 
D and reduce by unity the in-degrees of all its nodes in S, then reduce by 15*1 the 
out-degree of node i. The bds D' thus obtained will be called the bds reduced by S 
about node i from bds D. Equivalently, D' is the residual sequence obtained from D 
after connecting an out-stub from i to an in-stub of every node from S. 

Theorem 2 (Star-constrained graphicality) Let D be a bds in normal order on 
N nodes, and let Xi, \Xi\ < N — 1 — d ° , be a set of nodes whose in-stubs are forbidden 
to be connected to the out-stubs of node i (including i). Define Ci as the set of the 
first ("leftmost") d!f^ nodes in D but not from Xi. Then, there exists a simple digraph 
which realizes D and avoids connections from i to Xi, if and only if the bds D' reduced 
by Ci about node i from D is graphical. 

The proof of this theorem is found in Ref. [30]. What this theorem does is to turn 
a star-constrained graphicality problem for bds D into an unconstrained one on the 
reduced bds D'. The graphicality of D' is then easily tested via the FR theorem. The 
set Li as defined above will be called the leftmost set for node i. 

Although announced in its full generality, as Xi could be any predefined subset of 
nodes with \Xi\ < N—l—d^, this theorem applies directly to the digraph construction 
process when Xi represents the set of nodes to which connections were already made in 
previous steps from the same node z, hence forbidding us to make further connections 
from i to these very same nodes. In this case, the bds D represents the residual 
sequence D at that stage of the construction process. 

As discussed above, in order for us to be able to construct all the simple digraphs 
that realize a given bds, we need to find the allowed set Ai (D) for the next out- 
stub of i. Clearly, after every connection from the same node i, the residual sequence 
changes, and along with it the allowed set may change as well. In order to find Ai (D) 
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for the next out-stub of node i, we could just simply attempt connections sequentially 
to every node with non-zero in-degree not in Xi (D) U {i}, and test for graphically 
after each attempt using Th. [2] The set of nodes for which graphicality would have 
been preserved would form Ai[D^j . 

However, this would be inefficient and. actually, not needed. In fact, we can 
exploit a result which states that, if graphicality is broken by a connection, it will be 
broken by all other connections to the right of the previous one, in the normal order 
sense. This is expressed in the following: 

Theorem 3 Let D be a graphical bds in normal order and let Xi be a forbidden set 
for node i, with i G Afj. Let j < k be two nodes such that j, k £ Xi. Lf the residual bds 
Dj obtained from D after forming an edge directed from i to j is not graphical, then 
the bi-degree sequence D^, obtained from D by forming a directed edge from i to k is 
also not graphical. 

This theorem follows from the direct contraposition of Lemma 6 in Ref. [30] . Thus, 
what we need to do is to find efficiently the leftmost node q in the residual sequence 
in normal order, a connection to which would break graphicality. We will refer to this 
node q as the leftmost fail-node. All connections to this node and to nodes to its right 
are guaranteed to break (star-constrained) graphicality, whereas all connections to its 
left (with the exception of forbidden nodes and self) are guaranteed to preserve the 
graphical character. 

Note that both theorems [2] and [3] are based on the HH theorem for bi-degree 
sequences. In fact theorem [2] is a generalization of the HH theorem to include star 
constraints. Also note that, while for the FR theorem only the in-degrees must be 
ordered non-increasingly, for the HH theorem and hence for both theorems [2] and 
[3] the bds must be in normal order, as ordering by in-degrees only is not sufficient. 
This is easily seen from the following example of graphical bds (not in normal order) 
D 7 = {(2, 0), (2, 1), (0, 1), (0, 2)}. Using the HH theorem, if we do not worry about 
normal ordering, but just order by in-degree, we could choose to connect the out- 
stub of node (0,1) to an in-stub of node (2,0), then the out-stub of node (2,1) to 
the remaining in-stub of (2, 0) (connecting to the largest residual allowed residual in- 
degree), after which we have clearly broken graphicality: both out-stubs of (0, 2) now 
must be connected to the two in-stubs of (2, 1). 

We are now ready to present our digraph construction algorithm, which produces 
random samples from the set of all possible simple digraphs realizing a given bds. 

3. The algorithm 

Given a graphical bi-degree sequence D in normal order (initially D = D): 

1) Define as work-node the lowest-index node i with non-zero (residual) out-degree. 

2) Let Xi be the set of forbidden nodes for the work-node, which includes i. nodes 
with zero in-degrees and nodes to which connections were made from i. previously. 
In the beginning, Xi includes only the work-node and zero in-degree nodes. 

3) Find the set of nodes, Ai that can be connected to the work-node without breaking 
graphicality. 

4) Choose a node m G Ai uniformly at random and connect an out-stub of i to an 
in-stub of m. 

5) After this connection add node m to Xi. 
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6) If node i still has out-stubs, bring the residual sequence in normal order, then 
repeat the procedure from 3) until all out-stubs of the work node are connected 
away into edges. 

7) If there are other nodes left with out-stubs, reorder the residual degree sequence 
in normal order, and repeat from 1). 

The most involved step of the algorithm is finding the allowed set (step 3)), which is 
described next. 

3.1. Finding the allowed set 

Let i be the work-node chosen as in 1) and let D denote the normal ordered, residual 
sequence obtained after having connected some of the out-stubs of i to in-stubs of 
other nodes, such that graphicality is still preserved. These previous connections from 
node i form the set of forbidden nodes Xi for the next out-stub a of i. Xi also contains 
the work-node i itself i £ Xi and all other nodes with zero in-degrees. Let Ci be the 

set of the first (lowest index) d[°^ nodes from D, not in Xi. As D is (star constrained) 
graphical, we can connect a to any of the nodes in Ci without breaking graphicality 
(due to Theorem [2| , hence d C A4 . 

Let m be the last element of in the normal ordered bds D and let us "color" 
(label) red all the non-forbidden nodes, i.e., all the nodes not in A^, to the right of 
node m. Please note that these color labels are associated with the nodes, defined by 
their bi-degrees, and not with their indices of location in the sequence. This set of 
red nodes IZi forms the set of candidates for the leftmost fail-node q. All other nodes 
are colored (labelled) black. To find the leftmost fail-node we could simply connect 
out-stub a to an in-stub of a red node £, add the new connection temporarily to the 
set of forbidden nodes, bring the new residual sequence into normal order, then test 
for graphicality using Theorem [2] This procedure could then be repeated sequentially, 
with I going over all the red nodes from left to right, until graphicality would fail for 
the first time at i = q. However, the considerations in the following paragraphs allow 
us define a better method. 

For the sake of argument let us perform the sequential testing as explained above. 
It would imply the following steps for a given red node I : 

(a) Reduce the out-degree at the work-node i and the in-degree at I by unity, that is 
d° 1 — y dj — 1 and d% 1— >• d^ — 1, resulting in a new residual bds D^. 

(b) Bring into normal order (required by Theorem 2). Note that £ is the only node 
whose in-degree has changed and only the work-node had its out-degree changed 
(its in-degree was not affected). Thus, when bringing into normal order, the 
relative positioning of all the other nodes does not change. The work-node might 
have shifted to the right to a new position i' within the block of nodes with the 
same in-degree {if > i), and the red node's new position £' might have also moved 
to the right in the normal ordered sequence (£' > £). 

(c) Add I' to the forbidden set for the work-node. 

(d) Now, as required by Theorem [2] reduce by unity the in-degrees of the nodes in 
the left-most adjacency set Cy, and reduce the out-degree of the work-node i' to 
zero. This results in the new sequence DV- 

(e) Order the bds DV by in-degrees, non-increasingly. 

(f ) Apply the FR theorem to test for graphicality. 
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Thus, whether the connection of the work-node i to £ breaks graphically, ultimately 
depends on whether the residual bds DV fails (or passes) the FR test. However, as we 
noted before, for the FR test we do not need to have the bds DV in normal order, we 
only need to have it ordered non-increasingly by the in-degrees. Additionally, observe 
that in step (d) the reduction of the in-degrees always happens on the same set of 
nodes, independently of the red node £, that is the left-most adjacency set d> is the 
same for all £. Thus, in this particular case of Theorem 2's application, ultimately 
we do not need to bring into normal order (step (b)), only non-increasingly by 
in-degrees, which would be done anyway in step (e). That means we can just skip 
step (b), we do not need to move around any of the nodes at that stage. Thus, the 
only difference between the sequences DV for different £-s is at the position of this 
node after the reordering in (e), with respect to the rest of the sequence. 

These observations suggest that we should define a bds D' obtained from the 
bds D by reducing by unity the in-degrees of all nodes in the set Ci \ {m} and by 

d[°^ — 1 the out-degree of the work-node i, leaving only one out-stub (out-stub a) at i. 
Clearly, the bds D' is graphical (connecting out-stub a to an in-stub of node m surely 
preserves graphicality, by Theorem [2]) . Let us now order D' non-increasingly by its 
in-degrees, in a specific way, described as follows. Shift only the reduced in-degree 
nodes in D to the right with respect to the rest of the sequence such as to restore 
non-increasing ordering by the in-degrees (if needed) . Since only the in-degrees of the 
nodes in the set Ci \ {m} have been reduced, keep the relative ordering of all other 
nodes in D' exactly the same as in D. Thus the relative ordering of the red nodes and 
of the work node have been preserved as well. Let us denote the new location of the 
work node in D' by j (J < i) . Connecting now a to an in-stub of a red node £ in this 
sequence will produce the same set of residual bi-degrees as in step (d) above. To be 
able apply the FR theorem, then all we need to do is to shift to the right node £ in the 
sequence (if needed) to make sure that it is non-increasingly ordered by in-degrees. 
Since only the in-degree at £ was modified (reduced by unity) , this reordering is very 
simple: if x denotes the location of the last node of the block of nodes with the same 
in-degree as node £ in D' (x > £), then we simply swap the node at £ with the node 
at x after the reduction of the in-degree at £. Let us denote the obtained sequence by 
D". Clearly, it is non-increasingly ordered by in-degrees, and thus we can apply the 
FR theorem to see if it is graphical. Note: it could happen that x = j (e.g., there are 
many nodes with zero out-degree but larger in-degree than the work-node as defined 
in 1)), however, the steps below can be applied just the same. 

Next, we show how to identify the leftmost red fail-node q by investigating how 
the inequalities in (J3J) break down. Since D' is graphical, we have for all 1 < k < n — 1 
(n is the last element of D') that L'{k) < R'(k), where L' and R' are the left hand 
side (lhs) and the right hand side (rhs) of inequalities ^ written for D': 

L'(fc) = £df, (4) 

k n 

R'(k)= ^mmh-l^l+^mmh^Y (5) 

s=l s=k+l 

Let us denote by L"(k) and R"(k) the lhs and rhs of the inequality ^ corresponding 
to D". Since the rhs of ^ involves only out-degrees, and we only reduced the out- 
degree of the work-node from 1 to 0, we will always have R"(k) = R'(k) — 1, except 
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when k = 1 and the work-node is at j — 1, in which case R"(l) — -R'(l)- However, 
in this case, L"(l) = L'(l), because only the in-degree of j = 1 appears, which does 
not get changed. Thus, since L'(l) < R'(l) (D' is graphical), graphicality cannot be 
broken at fc = 1 when j = 1. Let us now consider that the work-node is still at position 
j = 1, but fc > 1. For 1 < k < x, the in-degrees in D" are the same as those in D', 
hence L"(k) = L'(k). For k > x, however, we have L"(k) = L'(k) — 1. Now consider 
j > 1. For 1 < fc < x, we have L"(fc) = L'(fc) and for fc > x, L"{k) =L'(k)-l. The 
following summarizes the relationships above: 



(A) j = 1: 



(A.l) fc = 1: 




= L'{1) , 


R"(l) 


= R'(1). 




(A.2) 1 < fc < x: 


i"(fc) 


= L'(k) , 


R"{k) 


= R'(k) - 


1 


(A.3) x < fc: 


L"(k) 


= L'(fc)-l, 


R"{k) 


= R'(k)- 


1 


J > 1: 












(B.l) l<k <x: 


L"{k) 


= L'(k) , 


R"{k) 


= R'(k) - 


1 


(B.2) a; < fc: 


L"(k) 


= L'(k) - 1 , 


R"(k) 


= R'(k) - 


1 



Since L'(k) < R'(k) for all fc, graphicality for D" can only be broken (that is to 
have L" > R" for some fc), if L'(k) = R'(k), namely in cases (A.2) and (B.l) above. 
Observe that L'(fc) and i?'(fc) are computed from D', hence they arc independent from 
£ or x. This gives us the following simple procedure for finding the leftmost fail-node, 
if it exists. Starting from k — 2 for j = 1, and fc = 1 for j > 1, find the smallest 
fc for which L'(k ) = R'(k ). If no such fc exists, then there are no fail-nodes and 
all non- for bidden nodes are to be included in the allowed set. If there is such a fco, 
the first red node q' to the right of fco (q' > fco + 1) is the leftmost fail-node of D", 
which when identified in the original bds D will give the leftmost fail-node q. All 
non- for bidden nodes to the left of q are to be included in the allowed set. 



3.2. Summary for finding the allowed set 

What we discussed in detail in the previous subsection corresponds to step (3) of the 
main algorithm described in the beginning of Section [3] Given the normal-ordered 
bds D at the end of step 2) of the main algorithm: 

(3.1) Identify Li from the first d{°^ nodes not in Afj. 

(3.2) Identify the "red" set IZi as those nodes that are neither in Li nor in A?,. Note, 
the color label is associated with the node, not its index. 

(3.3) Build D 7 as follows: 

= f df 3 - 1 ^ b G A \ {m} 
| d[ otherwise 

and 

_ (o) ( 1 if c = i 
c | d^ otherwise 

where m is the last node in 

(3.4) Shift nodes from Ci \ {m} to the right in the sequence (and only these) such as to 
restore ordering non-increasingly by in-degrees (if needed), preserving the color 
labels of the nodes in the process. The work-node may have shifted to a new 
location j after this step. This is the updated sequence D'. 
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(3.5) Starting from k = 1 if j ^ 1 or from k = 2 if j = 1, find ko as the smallest k 
such that L'(k) — R'(k), where L'(fc)and R'(k) are computed from the reordered 
(after step (3.4)) D' using ([i} and |H|. If there is no such ko, then the allowed 
set Ai is all the nodes in D except nodes from the forbidden set Xi. 

(3.6) Otherwise find the leftmost red node q' in the updated bds D' to the right of ko, 
that is with q' > fco . Then the corresponding node q in D , will be the leftmost fail 
node. Note that q' is the new position of the node at q in D after the reordering 
in (3.4). 

(3.7) The allowed set Ai is formed by all nodes in D not in Xi, and to the left of q. 



4. The sampling problem 

The algorithm generates an independent sample digraph every time it runs, without 
restarts or rejections, and it guarantees that every possible realization of a graphical 
bds by simple digraphs can be generated with a non-zero probability. However, the 
algorithm realizes the digraphs with non-uniform probability. Nevertheless, knowing 
the relative probability for every digraph's occurrence allows us to calculate network 
observable averages as if they were obtained from a uniform sampling. In particular, 
the following expression, which is a well-known result in biased sampling |43[ [44] . 
provides these averages as: 

V\^Li io(s,*)Q(s,) 

(Q) - j " m - (6) 

where Q is an observable measured from the samples Sj generated by an algorithm. 
The w(sj) sample weight is the inverse of the relative probability of the occurrence of 
Sj and M is the number of the samples generated. In Subsection 4.1 we give a detailed 



derivation of this formula, specialized to our graph construction problem. The weights 
of the samples generated by our algorithm are given by 

4 0) 

w(s) = Uf[k t (j). (7) 

i j=l 

where i runs over all the nodes with non-zero out-degree as they are picked by the 
algorithm to become work-nodes, and ki(j) — \Ai(j)\ is the size of the allowed sets 
Ai(j) J us t before connecting the j-th out-stub of i. Note that w > 1 since there always 
exists at least one digraph realizing the bds. Subsection 4.2 gives a derivation of d7l). 



4-1. Biased sampling over classes 

Our algorithm sequentially connects all stubs from a series of work nodes and finishes 
with a simple, labeled digraph. This process can be uniquely described by a path of 
connection sequences. Having chosen a work node i% for the first time, it determines 
the allowed set Ai x . We next choose uniformly at random a node ii(ii) € Ai 1 and 
connect a stub of i\ to a stub at Ji(ii). We could have chosen j\(i\) following any 
other criterion, but in that case the expression of the weights would have to be 
modified accordingly. After this connection we recompute the new allowed set A 3l (i\), 
then connect another stub of i\, and so on until all the stubs have been used up at i\. 
Let us denote by s such a path of connection sequences: 

s = I ii,ji{ii), ■ ■ ■ J^o){i 1 );i 2 ,j2{i2), ■ ■ ■ , j^ife) ■■■> (8) 
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where c^°' denotes the residual out-degree of node i. A path s uniquely defines the 
digraph G(s) created, as the collection of all connections in ^ forms the edge set of 
the created graph G(s). However, several paths may lead to the same digraph. Also 
note that the order of the connections in ^ matters in the calculation of the weight, 
as the corresponding allowed sets in general depend on history of connections up to 
that point. For a finite bi-degree sequence the number of distinct samples (paths) is 
also finite. Let us denote this set of paths by: 

II = {ttx, . . . ,7T P } . 

Let us now assume that we built with our algorithm a sequence of samples 
Si, s 2 , . . . , sju, and that the sample number M is large enough for us to see all elements 
of II sufficiently many times. Given some path s we compute a quantity Q(s), and 
we are interested in calculating the average of Q over path ensembles. In our case Q 
is defined on the final graph itself Q(s) = Q(G(s)), but for now we will not consider 
that, explicitly. If we were just simply computing the average of Q over the set of 
samples, we would obtain an average biased by the way the algorithm builds the paths 
from II: 

1 M P M 

W) = mE*) = E^*^ ( 9 ) 

1=1 fc=l 

where Mk is the number of times we have seen path tt^ appear in the sequence of 
samples. Clearly, 

p k = lim ^ (10) 

is the probability by which path Hk is generated via the algorithm. We now assume 
that we can compute analytically the path probabilities pk, from knowing how the 
algorithm works. Instead of (|9| we want to compute the average as if it was measured 
over the uniform ensemble of paths, that is: 

i p 

fc=l 



If we form: 



_ Ej=i p (bqQ( s <) , , 

\Q)bp - — M i — y 12 > 

2—ii=\ p(si) 



fc=l Mp(lT k ) 



we have limjvf->oo(Q)&}> = (Q)up, due to (10 1. Thus, the weighted average (12 1 should 



be used in order to obtain averages according to uniform sampling in the M 3> 1 limit . 

Let us assume that there is an equivalence relation "~" between paths, hence 
inducing a partitioning of II into K equivalence classes: II = C\ U . . . U Ck, where 
Ci = (ttj^, . . . ,7Tj.« \. The size of class Ci is denoted by pg — \Cg\. We have 

Y]f—i = P- Alternatively, for some given path tt, we will denote by C(ir) the 
equivalence class of tt and by h(tt) — \C(tt) \ its size. Let us also assume that if s, r £ Ci, 
that is s ~ r, then Q(s) — Q(r). For example, in our case distinct paths may lead to 
the same digraph. We introduce the equivalence relation "~" and say that two paths 
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s and r are equivalent, s ~ r if they produce the same labeled digraph, G(s) = G(r). 
Clearly, if Q depends only on the constructed graph, i.e., Q(ir) — Q (G(tt)) for all 
7r € II, then Q(s) — Q(r) whenever s ~ r. 

Our goal is to obtain the average of Q uniformly over the equivalence classes, that 

is: 



1 K 

(Q)uc = j?^Q i^ki) , 



(13) 



e=i 



where we chose to write the first element of Ci in the argument of Q, but of course, 
any other element could have been chosen from the same class, as Q is constant within 
a class. In general, (12) will not produce (<5) uc , but a sum weighted by class sizes. 



Instead, let us consider: 



(Q)bc = 

It is then easy to see that 

(Q)bc = 



M l 

1 fJ.(Si)p(Si) 



2^k=l M(7r fc )p(^)^ 1 fej Af- 



(14) 



Efc=i 



M k /M 



In order for (14 1 to be useful in practice, one has to be able to compute the size 



of the equivalence class /x(s) from seeing s and knowing how the algorithm works. 
Fortunately this is possible in our case, as shown next. 



4-. 2. Computing the weights 

First, let us note that when connecting the out-stubs of a work-node we are not 
affecting the out-stubs of any other nodes, but only in-stubs. Hence, all nodes with 
non-zero out-degrees will eventually be picked as work-nodes by the algorithm. Since 
normal ordering is first by in-degrees, the order in which nodes will become work- 
nodes depends on the sequence of connections. Let us now calculate the probability of 
the path s in Given a residual sequence, the work-node i\ is uniquely determined 
by the algorithm as described before. Since the next connection is picked uniformly 
at random, the probability of the link from i\ to j\(i\) £ is \Ai 1 (ii)| • Let 

ki{j) = \Ai(j)\ denote the number of nodes in Ai(j). Then, it is easy to see that the 
probability of a path s is given by: 



P(s) 



fc j=i 



(15) 



where i\, ■ ■ . , denote the work-nodes in the order in which they are picked by the 
algorithm. This expression can be computed readily in a computer as the algorithm 
progresses. In order for us to use (14 1 it seems that we would need also to obtain the 



size /x(s) of the class to which path s belongs. Clearly, two different paths s and s' 
will result in the same graph (s ~ s') if and only if the sequence of connections in 
one path is a permutation of the connections in the other path. Hence, the class size 
/i(s) is nothing but the number of permutations of the connections, which is the same 
for all paths, that is, all classes have the same size /i. Since all connections are made 
from a node first before moving on to another, we have \i — YliLi ■ However, 
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Figure 3. Biased sampling on the example bds Dg. The measure monitored is 
Newman's assortativity coefficient r 1451 . In b) the ensemble average was taken 
over 50 runs. 



we actually don't need to use this number: one can simply multiply by /i both the 
numerator and the denominator of (14 1 to obtain (6][7). 



4-.S. A simple example 

In this subsection we illustrate the algorithm on a simple sequence: Dg = 
{(2, 2), (2,1), (1,3), (1,1), (1,0)}. There are 11 distinct labeled digraphs realizing 
this sequence and there are 2!l!3!l!0! = 12 paths in a class, leading to the 
same graph. Two paths that lead to different graphs are for example si = 
{(1,4)(1,2);(3,1),(3,5),(3,2);(2, 1);(4,3)} (connect an out-stub of node 1 to an in- 
stubofnode4, etc.) and s 2 = {(1, 2)(1, 3); (2, 1); (4, 1); (3,4), (3, 5), (3, 2)}. For the 
former, iu(si) = [p(si)] _1 = 8 and for the latter it is w(&2) — [p^)] -1 = 54. Let us 
now consider the Pearson coefficient r of degree-degree correlations, or the assortativity 
coefficient defined for directed graphs [IS] as our network observable Q — r. For each 
one of the 11 graphical realizations of Dg,r can be calculated exactly, as can the 
uniform average over this ensemble, obtaining {f)theo — 

-0.040506. We will refer to 

(r)theo as the "theoretical value". We then let our algorithm run on this sequence 
to produce M samples and using (6][7) to obtain the corresponding coefficient (r)jvf- 



Fig |3p) shows a few runs with different seeds and their convergence to the theoretical 
value. Fig^) shows the standard deviation {[(r)M ~ (r)theo] 2 ) 1 ^ 2 where the overline 
denotes an ensemble average over runs. 

5. Complexity of the algorithm 

To determine the theoretical upper bound for the complexity of the algorithm, that 
is the worst-case complexity, notice that there are only three steps in the algorithm 
that require more than O (1) computational operations, or steps, to complete. 

First, after each connection is placed, one must bring the residual sequence into 
normal order, steps 6) or 7). To accomplish this, both the work-node i and the target 
node m will have to move to the right, but the relative positions of all other nodes 
will remain unchanged. In other words, if we were to remove nodes i and m, the 
rest of the bds would already be sorted. Thus, in order to complete these steps, one 
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only has to find the new positions of nodes i and m and insert them into the already 
sorted bds. Therefore, the complexity of either one of step 6) and step 7) is simply 
O (21ogiV + N) ~ O (N), where N is the number of the nodes in the sequence being 
ordered. 

Second, the allowed set A must be built before placing each connection (step 3). 



Following the summary of this step, given in Subsection 3.2 notice that steps 3.1 
to 3.4 can be all finished during a single scan of the residual bds. This is clearly so 
for the creation of the leftmost set Li and for setting the "red" color labels (or flags) 
(steps (3.1) and (3.4)). Concerning the ordering of the bds D', it is possible to create 
it already sorted by simply scanning the bds D while keeping track of the in-degree 
d* of the nodes currently being copied and the index a in D' of the first node with 
that in-degree. Then, because D is in normal order, the only possibility for a node in 
D' to break the order is if its in-degree equals d* + 1. In this case, it can be simply 
swapped with the node at a, because, as argued in Subsection |3.1[ the mechanism to 
build the allowed set is entirely based on the FR theorem, which does not require the 
bds to be in normal order, but to be simply ordered non-increasingly by its in-degrees. 
Thus, steps (3.1) to (3.4) can be completed in O (N) steps. 

Third, the computation of the sums L' , R' and their comparison must be 
conducted, which is the same step as Q in an FR test. To determine the complexity of 
an FR test note that computing the repeated sums for each one of the inequalities (J3j) 
is quite inefficient. Instead, below we derive recurrence relations that allow us to 
complete the FR test in a linear, O (N) number of steps. 

The steps of the main algorithm are done sequentially, and thus can all be 
completed in a total of O (N) steps. They must, however, be repeated for each edge 
in the digraph. Thus, the maximum complexity of the algorithm is O (NM) where 
M = d\°^ is the number of edges. Since O (M) < O (N 2 ) the maximum complexity 
of the algorithm is 0(N 3 ). It is important to note though, that for a given bds the 
complexity of the algorithm can be substantially smaller, similar to the case for our 
undirected graph sampling algorithm 

5.1. The Fulkerson-Ryser test revisited 

The most complex part of the Fulkerson-Ryser test is to compute the lhs and the rhs 
of inequalities pi, which we rewrite here for the sake of readability: 

k 



s=l 



A' 



R(k) = ^minjfc- l,di o) }+^min{fc,4°)}. 

s=l s=k+l 

Our goal is to find recursion relations for L(k) and R(k). For the lhs the relation 



is simply 



L(k + 1) = L{k)+d { * } 



with L(l) = d[ l) . 

For the rhs, first note that one can write it as 

N 

R{k) = -/c + ^min{fc,. 9j (o) (fc)} , (16) 



i=l 
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where gf\k) is the family of integer sequences defined as 

£\k)- 



d[ o) + 1 Mi C fc 



d[ o) Mi>k 

Now, let us introduce G k (p) = $Zi=i ^ p 9 <o) (fc)' that is, the number of indices i for 
which g\°\k) = p. Then, from {l6| follows that 

k N 
R(k) = -k + ^ P G k (p) + k J2 G * (p) , (17) 

p—1 p=k+l 

hence 

R(l) = N- 1-Gi (0) , (18) 

where we used the fact that X^^Lo Gk (?) = ^ ■ 

Furthermore, let us introduce the following notations: 

AG fe ( P ) = G k ( P ) - G fe _ x (p) 
G k {q) = Y^G k {i) . 



Then, after some simple manipulations, from (17) it follows that 
R(k) -R(k-1)=N-1- G fc _i (fc - 1) 

fc-1 N 

+ ]>>AG fe (p) + fc^AG fc (p) . (19) 

p=l p=k 



Finally, notice that AG k Op) = 5 ,< ) . — 5 ,(<,). Substituting it into (19 1, we obtain: 

R(k - 1) + N - G k _!{k - 1) Vd k o) <k 

/,'(/,) = <; (20) 

R(k - 1) + N - G k -i(k - 1) - 1 Vd k o) ^k 

Thus, we have turned the problem of finding a recursion relation for R(k) into 
the problem of finding G k (k). To solve this, first note that 

G k (k) = G fe _! (k - 1) + G fe _! (fc) - <5 M <o, , 

with Gi (1) = Gi (0) +Gi (1). The above equation constitutes a recursion relation for 
G k (q). Such a relation can be rewritten as 

G k (fc) = G fe _! (k - 1) + d (A;) + S (fc) , 

where 

fc-i fc 

t=2 {=2 

Observe that S (fc) and G\ (fc) can be easily computed while scanning the bds, and 
then calculating L(k) and R(k) for each fc requires a single operation. Thus, the entire 
FR test can be completed in O (N) steps. 
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300 



log w 



Figure 4. Probability distribution p of the logarithm of weights for an ensemble 
of bi-degree sequences on N = 100 nodes. The in-degrees were drawn from a 
normalized power-law distribution ~ with 7 = 3 and the out-degrees were 
drawn from a Poisson distribution e~^\ d ° ut /d out \, with the same average as the 
average in-degree, A = (din)' The black circles are the simulation data and the 
red continuous line is a Gaussian fit. 



6. Discussion 

In summary, we have developed a graph construction and sampling algorithm to 
construct simple directed graphs realizing a given sequence of in- and out-degrees. 
Such constructions are needed in practical modeling situations, ranging from epidemics 
and sociology through food-webs to transcriptional regulatory networks, where we are 
interested in learning about the statistical properties of the network observables as 
determined only by the bi-degree sequence and nothing else. 

Unlike existing algorithms such as the Configuration Model, which is affected by 
uncontrolled biases and unacceptably long running times except for a very restricted 
class of sequences, our algorithm is rejection-free. Also, it guarantees the independence 
of the produced samples, unlike MCMC methods, which have unknown mixing 
times. While its mathematical underpinnings are nontrivial, the algorithm itself is 
straightforward to implement. In principle, our approach can be extended to include 
more complex constraints, such as a given sequence of motifs frequencies, but we have 
only concentrated on degree sequences since they are, arguably, the most fundamental 
of constraints. The algorithm can also be used to sample from given in- and out- 
degree distributions, not just sequences: given such distributions, one first samples a 
graphical bds from these, then one applies our algorithm to generate digraphs. In this 
case, however, the sample weights ^ must be modified to reflect the probability of 
the occurrence of the given graphical bds when drawn from the distributions. 

Just as in the case of undirected graphs, we can expect the distributions of the 
weights for large graphs to be log-normal, as shown in Ref. [11J. As an example, 
figure [4] shows the distribution for bi-degree sequences in which the in-degrees follow 
a power law with exponent 7 = 3 and the out-degrees a Poisson distribution whose 
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Figure 5. Mean m (black circles) and standard deviation a (red squares) of 
the distributions of the logarithm of the weights vs. number of nodes N of 
samples. In-degrees and out-degrees are both drawn from a power-law distribution 
P (d) ~ d~ J , with 7 = 3. The solid black line and the dashed red line are data fit 
results, showing that m and a follow power-law scaling laws m ~ TV" and a ~ TV 3 . 
The values of the exponents, given by the slopes of the lines are a = 1.23 ± 0.02 
and $ = 0.81 ± 0.02. 

mean matches the average in-degree. Indeed, the distribution of the weight logarithms 
is well approximated by a Gaussian. Similarly the undirected case, we find for all the 
examples we studied numerically, that the standard deviation a of the distributions 
of weight logarithms grows slower than the mean m with the number of nodes N; 
see figure [5] showing the scaling of m and a for bi-degree sequences in which both in- 
degrees and out-degrees follow a power law distribution with exponent 7 = 3. Thus, 
we may expect that typically, in the N — > oo limit, the rescaled weight distribution 
becomes a delta function, making the sampling asymptotically uniform. 

Bounds on the complexity of the algorithm could easily be obtained by inspecting 
the algorithm, showing a maximum complexity on the order of O(NM) where M is 

the total number of edges, M = Y^Li <H • 

In developing our results, we also provided an efficient way of implementing the 
Fulkerson-Ryser test, whose scope of application goes beyond our present algorithm, 
as it can be used in any context to determine whether a bi-degree sequence is graphical. 
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